Linear regression for calibration curve fitting

GCCL에서 제공되는 분석서비스는 정말 다양합니다. Small molecule 분석에 사용되는 질량분석을 비롯해서 Large molecule 분석에 사용되는 LBA, PCR은 물론, 여기에 다양한 진단 검사까지 포함하면 일일이 열거하기 힘들 정도로 많습니다. 그러나 이런 다양한 분석도 대부분은 정량 분석(quantitative analysis)이 목적입니다. 정량 분석하면 Calibration curve fitting을 빼놓을 수 없습니다. 분석 업무에 종사하거나 분석 결과를 다뤄야 하는 업무중이시라면 늘 주의깊게 살펴보는 작업일 것입니다. 그런데 curve fitting에서 사용되는 수학적 원리를 설명하라면 그런게 있었어? 라며 갸우뚱 하실 것 같습니다. 우리는 보통 calibration curve를 fitting한다는 표현을 하는데 구분해 보면 Calibration curve가 직선의 형태일 때 일반적으로는 Linear regression curve fitting이라고 합니다. 직선이 아닌 경우는 Non-linear regression curve fitting이라고 표현합니다. 여기서 regression은 회귀라는 뜻입니다. 따라서 curve fitting은 곧 회귀 분석입니다. 또는, “왜 임상시험검체 분석을 위해 분석법을 개발하고 평가할 때 Calibration standard curve의 $R^2$ 은 평가하지 않나요?”
라는 질문을 받기도 합니다.
우리가 많이 보는 규정에도 이를 평가하라는 말은 없습니다. 왜일까요?
사실 평가 할 필요가 없어서 그런 것인데 그렇다면 왜 필요하지 않을까요?
답을 하기 위해서는

R^{2}

$R^2$ 의 의미를 정확히 알아야 하는데 바로 이런 부분에서 수학과 통계학이 필요합니다.
따라서 필요한 수학적 또는 통계학적 개념과 원리를 잘 알아둘 필요가 있습니다.

이런 이유로 앞으로 몇 차례에 걸쳐 수학적 또는 통계학적 개념 중에서
임상시험검체 분석에 사용되는 몇가지 중요한 개념을 심도있게 다뤄보고자 합니다.
오늘은 그 첫번째로 선형 회귀에 관한 내용을 다뤄보겠습니다.
계산과 설명을 위해 R이라는 통계소트웨어가 사용되었다는 점을 알려드립니다.

선형 회귀(Linear Regression) 분석은 통계학과 머신러닝에서 널리 사용되는 기본적인 방법입니다.
독립변수(예측변수)가 종속변수(반응변수)에 미치는 영향을 모델링하는 데 사용하는데 이름에서도 알 수 있듯이
선형 회귀는 이 두 변수이 관계가 직선 형태로 나타난다는 전제로 수행됩니다.

1.1 선형 회귀 일반

선형 회귀 방정식은 아래와 같은 형태입니다.

y = b_{0} + b_{1} x + ϵ

$y = b_0 + b_1 x + \epsilon$

각 기호가 의미하는 바는 이렇습니다.

$y$ : 종속변수(예측하려는 값. 결과변수라고도 불림)
$x$ : 독립변수(예측에 사용되는 값. 예측변수라고도 불림 )
$b_0$ : 절편(intercept)
$b_1$ :기울기(slope)
$\epsilon$ : 오차 항

위와 같은 수식 자체를 회귀 모형(Regression model)이라고 부르는데 회귀 분석은 주어진 데이터에 맞는
최적의 회귀 계수 (Coefficients 또는 Parameter라고 부르고 여기서는 $b_0$ , $b_1$ 을 의미)를 찾는 행위라고 할 수 있습니다.

1.1 선형 회귀의 종류

선형 회귀는 크게 두 가지로 나눌 수 있습니다.

단순 선형 회귀 (Simple Linear Regression): 하나의 독립변수를 사용하여 하나의 종속변수를 예측할 때 사용합니다.
다중 선형 회귀 (Multiple Linear Regression): 여러 개의 독립변수를 사용하여 하나의 종속변수를 예측할 때 사용합니다.

예측한다는 표현이 낯설 수 있는데 회귀 분석은 곧 예측 모형을 만드는 작업이므로 이렇게 예측한다는 표현이 자주 등장합니다.

1.3 선형 회귀의 일반적인 가정

선형 회귀에 필요한 몇 가지 기본적인 가정들을 요약하면 다음과 같습니다.

선형성 (Linearity):
변수 간의 관계가 선형적이어야 합니다.
이 뜻은 독립변수의 값이 증가함에 따라 종속변수의 평균값들이 하나의 직선을 형성한다는 뜻입니다.
아래 그림은 linear calibration curve fitting을 시뮬레이션 해본 것입니다.
각 개별 분석에서 얻어진 데이터들을 모아서 표현하면 아래 그럼처럼 보이게 되는데 Calibration curve (여기서는 직선이지만)가
각 calibrator 수준에서 데이터의 평균을 지나고 있는 것을 볼 수 있습니다.
종속변수의 평균값들이 하나의 직선을 형성한다는 뜻은 이러한 개념이라고 생각하시면 됩니다.

오차의 독립성 (Independence):
관측치 간에 독립성은 임의의 두 관측값의 잔차(Residual)들이 무관(서로 독립적)해야 한다는 뜻입니다.
다른 말로는 자기상관(Autocorrelation)이 없다고 표현합니다.
여기서 말하는 잔차는 모형이 예측한 값과 실측값의 차이를 의미합니다.
등분산성 (Homoscedasticity):
독립변수들의 각 수준에서 잔차들의 분산이 일정해야 한다는 뜻입니다.
잔차는 작을 수록 좋은데 크기 뿐만 아니라 0을 기준을 모든 범위에서 고르게 분포하는 것이 좋습니다.
그 이유는 예측값과 실제 측정값이 똑같을 때 잔차가 0이기 때문입니다.
아래 그림은 위 Calibration standard sample의 분석 결과 중 잔차의 분포를 표현한 것입니다.
각 수준에서고르게 분포하는 모습을 볼 수 있는데 이때 잔차의 분포를 등분산이라고 표현합니다.

model <- lm(y ~ x, d)
resd <- resid(model)
plot(fitted(model), resd, main = "Residual plot",
     xlab = "Fitted values",
     ylab = "Residuals")
abline(0,0)

정규성 (Normality):
잔차의 분포가 정규 분포를 따라야 합니다.
이 가정은 모형의 잔차들이 평균이 0인 정규분포를 따르는 무작위 값들이라는 뜻입니다.
쉽게 말해서 회귀 모형과 관측값의 차이들이 대부분 0이거나 0에 가까운 값이며 크기가 0보다 훨씬 큰 잔차들은 가끔 나타난다는 뜻입니다.
아래 Density plot으로 보면 가운데(평균)가 0이며 잔차의 분포 자체가 가운데 0을 중심으로 밀집되어 있는 형태를 볼 수 있고
양 끝으로 갈 수록 그 출현빈도가 낮아지는 모습을 볼 수 있습니다.
Destiny plot이 아닌 Histogram으로 표현하기도 합니다.

plot(density(resd), main = "Density plot",
     xlab = "Residual")

hist(resd, main = "Histogram of Residual")

1.4 선형 회귀의 장단점

선형 회귀의 장단점을 살펴보면 이렇습니다.

장점:

이해와 구현이 쉽습니다.
계산이 효율적이며 빠릅니다.
해석이 직관적입니다.

단점:

선형성 가정이 성립하지 않으면 성능이 떨어집니다.
이상치(outlier)에 민감합니다.
다중공선성(Multicollinearity)이 있을 경우 문제가 발생할 수 있습니다.

다중공선성은 회귀 모형에 여러 개의 예측변수가 사용될 때 한 예측변수가 다른 예측변수와 높은 상관관계가 있는 경우를 의미합니다.
이 경우 두 예측변수의 고유한 회귀 계수 추정이 불가능하다는 문제가 생기는데 calibration curve fitting을 위한 선형 회귀에서는 발생하지 않는 문제입니다.

본격적으로 선형 회귀 분석의 대표적인 예인 linear regression curve fitting을 살펴보겠습니다.
Calibration curve fitting 결과가 $Y = Ax+B$ 의 형태가 될때 선형 회귀 분석이 사용되는데 이때
변수는 standard sample의 농도와 그 assay signal입니다. 따라서 물질의 농도는 예측변수 또는 독립변수가
되고 assay signal은 반응변수 또는 종속변수가 됩니다.
그러면 이때 예측변수의 숫자는 몇개일까요?
Calibration standard sample의 level(농도)만큼 많을까요?
예측변수의 종류는 농도변수 뿐이므로 예측변수는 1개 뿐입니다.

아래와 같은 Calibration standard sample의 측정 데이터가 있다고 가정해보겠습니다.
X는 농도, Y는 assay signal입니다.

df_data <- data.frame(x = c(1, 2, 4, 8, 16, 32, 64, 128), 
                      y = c(0.5, 1.1, 1.9, 4.0, 7.8, 16.1, 31.4, 68.1 ))

ggplot(df_data, aes(x, y))+
  geom_point()+
  theme_bw()+
  xlab("Concentration (ng/mL)") +
  ylab("Assay response")

이 데이터로 R을 이용해서 회귀 모형을 만들어 보겠습니다.

m <- lm(y ~ x, df_data)

만들어진 회귀 모형을 호출하면 아래처럼 보입니다.

summary.lm(m)

## 
## Call:
## lm(formula = y ~ x, data = df_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9099 -0.2229  0.2377  0.4435  1.0273 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.452982   0.422806  -1.071    0.325    
## x            0.527545   0.008091  65.200 8.75e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9476 on 6 degrees of freedom
## Multiple R-squared:  0.9986, Adjusted R-squared:  0.9984 
## F-statistic:  4251 on 1 and 6 DF,  p-value: 8.754e-10

화면 중간 Coefficients 중 Estimate가 우리가 필요로 하는 가장 중요한 정보인 Calibration curve의 기울기와 절편의 정보를 담고 있습니다.
그 정보를 이용해서 회귀 모형을 표현하면 아래와 같습니다.

y = - 0.4529822 + 0.5275445 x

$y= -0.4529822 + 0.5275445x$ Calibration curve를 표현하는 수식과 다르지 않습니다.
그러면 회귀 모형이 잘 만들어 졌는지 눈으로 확인해볼 차례입니다.

## regression curve
ggplot(df_data, aes(x, y))+
  geom_point()+
  stat_function(fun = function(x, y)(y = m$coefficients[2]*x-m$coefficients[1]), col="red")+
  theme_bw()+
  xlab("ng/mL") +
  ylab("Assay Response")

빨간 실선이 곧 회귀 모형이 예측한 값을 직선으로 연결한 형태가 되고 곧 Calibration curve가 됩니다.
여기서는 R을 사용했다는 차이만 있을 뿐, 우리가 사용하고 있는 분석소프웨어들 역시 비슷한 기능을 수행합니다.

그렇다면 어떻게 R을 포함한 소프트웨어들이 기울기와 절편을 계산할까요?
그 방법을 살펴보기 전에 한 가지 더 확인할 것이 있습니다.
회귀 모형의 목적은 예측이므로 예측 정확도가 높아야 좋은 모형입니다.
아래 데이터에서 y_hat 컬럼이 회귀 모형을 통해 예측된 fitted value(모형이 예측한 값)입니다.

df_data$y_hat <- m$fitted.values ## extraction of fitted values
df_data

##     x    y       y_hat
## 1   1  0.5  0.07456235
## 2   2  1.1  0.60210689
## 3   4  1.9  1.65719597
## 4   8  4.0  3.76737413
## 5  16  7.8  7.98773044
## 6  32 16.1 16.42844307
## 7  64 31.4 33.30986832
## 8 128 68.1 67.07271882

모형의 정확도를 판별하는 방법은 아래와 같습니다.

실측값과 예측값 간의 차이인 잔차(Residual)를 계산한다.
잔차는 부호(양수 또는 음수)가 존재하므로 부호를 상쇄하기 위해 제곱한다.
제곱한 잔차를 모두 더한 잔차제곱합을 계산한다.

제곱합은 각 잔차의 제곱을 합한다 뜻으로 잔차제곱합(Sum of squared residual)으로 부르기도 합니다.
약어로는 SSR 또는 RSS 또는 SSE로 표현하기도 하는데 여기서는 SSE로 사용하겠습니다.
수식으로 표현하면 아래와 같습니다.

S S E = \sum_{i = 1}^{n} (y_{i} - \hat{y_{i}})^{2}

$SSE = \sum_{i=1}^{n}(y_{i}-\hat{y_i})^2$ 앞에서 만든 회귀 모형의 SSE를 계산해보겠습니다.
위 수식으로 계산한 결과:

sum((df_data$y-df_data$y_hat)^2) ## manual computation of SSE

## [1] 5.387985

R이 계산한 결과:

sum(m$residuals^2) ## manual computation of SSE

## [1] 5.387985

R이 계산한 결과와 직접 계산한 결과는 같았는데 이론적으로는 SSE가 0인 모형이 가장 좋은 모형입니다.
SSE가 0이라면 실측값과 예측값의 차이가 전혀 없다는 뜻이 되므로 그 차이가 커지면 그 만큼 정확도가 떨어진다는 뜻이 됩니다.

3.1 Cost function

R을 비롯한 소프트웨어들이 선형 회귀 분석을 수행하거나 이를 이용해서 Calibration curve fitting을 할 때 사용하는 방법은
Ordinary Least squares method, 줄여서 OLS라고 부르고 우리 말로는 최소 제곱법이라고 합니다.

최소 제곱법을 살펴 보기 전에 한가지 알아야 할 개념이 있습니다.
Cost function(비용함수)에 관한 개념입니다.

아래 데이터는 기울기 2, 절편은 0인 단순 선형 회귀 모형의 데이터입니다.

x0 <- seq(from = 0 , to = 1, 0.1)
z <- 2*x0
df_cost <- data.frame(x0, z)
ggplot(df_cost, aes(x0, z)) +
  geom_point()+
  theme_bw()

모형의 예측값을 나타내는 수식은 아래와 같습니다.

\hat{y_{i}} = 2 x_{i}

$\hat{y_i} =2x_i$

기울기 2를 일반화해서 표현하면 아래와 같은 형태로 다시 쓸 수 있습니다.

\hat{y_{i}} = ω x_{i}

$\hat{y_i} =\omega x_i$

여기서는 절편이 0라서 회귀 모형의 정확도는 기울기( $\omega$ )를 어떻게 추정하는지에 따라 달라집니다.
그런데 여기서는 SSE가 아닌 mean squared error(MSE)라는 개념을 사용해보겠습니다.
아래 수식이 MSE를 표현하는 수식인데 잔차제곱합의 평균정도라고 이해할 수 있겠습니다.

J (w) = \frac{\sum_{i = 1}^{m} (y_{i} - w x_{i})^{2}}{2 m}

$J(w) =\frac{\sum_{i=1}^{m}(y_i-wx_i)^2}{2m}$

기울기를 0부터 4까지 0.1씩 변화시키면서 MSE가 어떻게 변하는지 살펴보겠습니다.

w <- seq(from = 0, to= 4, .1)
jw <- c()
for(i in 1:41){
  jw.0 <- sum(df_cost$z-w[i]*df_cost$x0)^2/(2*11)
  jw <- c(jw, jw.0)
}

df_jw <- data.frame(w, jw)
ggplot(df_jw, aes(w, jw))+
  geom_point()+
  ylab  ("J(w)")+
  theme_bw()

예상대로 $\omega$ = 2에서 J( $\omega$ )값이 0인 것을 확인할 수 있습니다.
위에서 MSE를 계산하는데 사용한 수식을 Cost function(비용함수) 또는 lost function(손실 함수)라고 부릅니다.
그리고 최소 제곱법은 비용함수가 0 또는 최소가 되는 회귀 계수를 계산하는 방법을 의미합니다.
즉 기울기에 2를 대입하면 SSE와 MSE가 최소가 됩니다.
주어진 데이터로 회귀 분석을 수행하면 이 사실을 확인할 수 있습니다(아래 참고).

summary.lm(lm(z ~ x0, df_cost))

## Warning in summary.lm(lm(z ~ x0, df_cost)): essentially perfect fit: summary
## may be unreliable

## 
## Call:
## lm(formula = z ~ x0, data = df_cost)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##      0      0      0      0      0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0          0     NaN      NaN    
## x0                 2          0     Inf   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0 on 9 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic:   Inf on 1 and 9 DF,  p-value: < 2.2e-16

따라서 회귀 분석을 통해 얻은 선형 회귀 계수는 바로 비용함수가 최소가 되는 조건을 의미합니다.
또는 최적화된 회귀 계수는 바로 MSE가 최소인 회귀 계수라고 할 수 있습니다.

그런데 우리가 처음 만든 회귀 모형은 기울기 뿐 아니라 절편도 포함하고 있습니다.
이 경우 비용함수는 3차원으로 표현되고 평면의 두 축인 x, y 각 축이 기울기와 절편, z축이 MSE가 되는 형태가 됩니다.
아래 interactive 3d plot이 바로 만들어진 모형의 cost function을 시각화한 결과물입니다.

Plot에서도 확인이 가능하듯이 비용함수가 최소가 되는 지점은 slope ~ 0.5, intercept ~ -0.5 지점입니다.
우리가 알고 있는 각 수치와 정확하게 일치하지는 않는데 기울기와 절편값을 0.1 씩 변화시켰기 때문입니다.
실제로 우리가 만든 모형의 비용함수의 최소값은 아래와 같습니다.

1/(2*8)*sum((df_data$y-m$coefficients[[1]]-m$coefficients[[2]]*df_data$x)^2)

## [1] 0.336749

3.2 Gradient Descendant method

일반적으로 비용함수를 최소화하는 선형 회귀 모형의 기울기와 절편은 최소 제곱법 또는 계산하거나 또는
Gradient descendant method(경사 강하법)으로 계산합니다.
경사 강하법은 최적화 알고리즘의 하나로 머신러닝에서 많이 쓰이는 개념인데 여기서 그 개념을 자세히 다 설명하기는 어렵습니다.
그러나 최소 제곱법과 구분되는 몇가지 특징이 있기 때문에 그 특성을 살펴보기 위해 아래에 해당 알고리즘을 구현했습니다.

gd <- function(x, y, start, max_iter, mu) {
  n <- length(y)
  alpha <- mu
  cost.his <- numeric(max_iter) 
  
  ls.b1 <- sum((x-mean(x))*(y-mean(y)))/sum((x-mean(x))^2)  
  ls.b0 <- mean(y)-mean(x)*ls.b1
  theta.ls <- c(ls.b0, ls.b1)
  
  
  cost.function <- function(x, y, start){
     m <- length(y)
     yhat <- x*start[2]+start[1]
     cost <- (1/(2*m))*sum((yhat-y)^2)
     return(cost)
     }
  
  for(i in 1:max_iter){
     p.b0 <- -(1/n)*sum(y-start[1]-start[2]*x)
     p.b1 <- -(1/n)*sum(x*(y-start[1]-start[2]*x))
     
     newb0 <- start[1]-alpha*p.b0
     newb1 <- start[2]-alpha*p.b1
     start <- c(newb0, newb1)
     
     cost.his[i] <- cost.function(x, y, start)
    
     
     if (all(abs(start - theta.ls) < 1e-6)) { 
      cat("Converged at iteration", i, "\n")
       cat(paste("B0:", start[1], "\n"))
       cat(paste("B1:", start[2]))
      break 
    }
    
    if (i == max_iter) {
      cat("Reached maximum iterations:", i, "\n")
      cat(paste("B0:", start[1], "\n"))
      cat(paste("B1:", start[2]))
    }
  }
  return(list(ls_theta = theta.ls, optimal_theta = start, cost_history = cost.his, iterations = i))
}

위 코드를 실행하면 최소제곱법과 경사 강하법으로 추정한 회귀 계수들를 비교할 수 있습니다.
먼저 절편과 기울기 추정에 사용된 초기값은 모두 0으로 보통은 Start value라고 표현합니다.

또 이 값들을 변화시킬 때 사용하는 parameter인 learning rate는 1e-6(0.000001)로 설정했습니다.
반복 횟수는 1000으로 제한하여 그 안에 우리가 원하는 답을 찾을 수 있는지 시험해봤습니다.

아래가 그 결과입니다.

start <- c(0,0)
x <- df_data$x
y <- df_data$y
result1 <- gd(x, y, start, 1000, 1e-6)

## Reached maximum iterations: 1000 
## B0: 0.005413044676279 
## B1: 0.488292087977798

아래는 최소제곱법으로 계산한 결과입니다. 앞이 절편, 뒤가 기울기입니다.

result1$ls_theta

## [1] -0.4529822  0.5275445

두 방법으로 계산한 절편과 기울기값이 다릅니다.
반복 횟수를 30만번까지 늘려보겠습니다.

result2 <- gd(x, y, start, 300000, 1e-6)

## Reached maximum iterations: 300000 
## B0: -0.0727155393416777 
## B1: 0.523104608022892

근접해가고 있는 것 같지만 아직 우리가 알고 있는 수치와는 거리가 있는 편입니다.
마지막으로 Learning rate를 수정해보겠습니다.

result3 <- gd(x, y, start, 300000, 1e-4)

## Converged at iteration 207644 
## B0: -0.452981184354306 
## B1: 0.52754452744114

드디어 알고있는 값을 찾았습니다. 이때 MSE도 위에서 계산한 MSE와 같습니다.

result3$cost_history[[result3$iterations]]

## [1] 0.336749

Start value 설정에 따라 다르긴 하지만 207644번의 반복과 learning rate의 조절이 필요했습니다.
따라서 linear regression curve fitting에 사용하기에는 매우 불편하다는 사실을 알 수 있습니다.

3.3 Least squares method

이제 최소 제곱법으로 기울기와 절편을 계산하는 방법을 살펴보겠습니다.
정규방정식(Normal equation)과 관련된 내용으로 조금 긴 유도 과정이 필요합니다.
참고로 그 유도과정은 아래 정리했으니 궁금하신 분들은 아래 왼쪽 화살표를 클릭 하시면 됩니다.

$b_1$ 에 관한 수식을 유도하는 과정

SSE는 아래와 같은 형태로 다시 정의할 수 있습니다.

S S E = \sum_{i = 1}^{n} (y_{i} - b_{0} - b_{i} x_{i})^{2}

$SSE= \sum_{i=1}^n(y_i-b_0-b_ix_i)^2$

SSE가 최소가 되는 회귀 계수는 결국 SSE가 0인 조건에 수렴하므로 위 식은 아래처럼 표현할 수 있습니다.

\sum_{i = 1}^{n} (y_{i} - b_{0} - b_{i} x_{i})^{2} = 0

$\sum_{i=1}^n(y_i-b_0-b_ix_i)^2=0$

이제 이 수식을 각 회귀 계수로 편미분합니다.
먼저 $b_0$ 으로 편미분하면

\frac{\partial S S E}{\partial b_{0}} = 2 \sum_{i = 1}^{n} (y_{i} - b_{0} - b_{i} x_{i})

$\frac{\partial SSE}{\partial b_0} = 2\sum_{i=1}^n(y_i-b_0-b_ix_i)$ 아래와 같은 수식을 얻을 수 있습니다.

2 \sum_{i = 1}^{n} (y_{i} - b_{0} - b_{1} x_{i}) = 0

$2\sum_{i=1}^n(y_i-b_0-b_1x_i) = 0$

$b_1$ 에 대해 똑같이 편미분하면

\frac{\partial S S E}{\partial b_{1}} = 2 \sum_{i = 1}^{n} (y_{i} - b_{0} - b_{i} x_{i}) x_{i}

$\frac{\partial SSE}{\partial b_1} = 2\sum_{i=1}^n(y_i-b_0-b_ix_i)x_i$

아래와 같은 수식을 얻을 수 있습니다.

2 \sum_{i = 1}^{n} (y_{i} - b_{0} - b_{1} x_{i}) x_{i} = 0

$2\sum_{i=1}^n(y_i-b_0-b_1x_i)x_i =0$

이렇게 얻은 두 식으로 연립방정식의 해를 구하는 방법을 이용해서 $b_0$ 과 $b_1$ 을 계산하는 수식을 다시 정의합니다.

먼저 $b_1$ 을 구하는 수식의 유도 과정을 요약하면 아래와 같습니다.

파란색 수식을 2로 나누고 양변에 $\sum x_i$ 를 곱하면 아래와 같습니다.

\sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i} - n b_{0} \sum_{i = 1}^{n} x_{i} - b_{1} (\sum_{i = 1}^{n} x_{i})^{2} = 0

$\sum_{i=1}^n x_i \sum _{i=1}^ny_i -nb_0\sum_{i=1}^n x_i - b_1(\sum_{i=1}^n x_i)^2 =0$

초록색 수식도 2로 나누고 양변에 n을 곱하면 아래와 같습니다.

n \sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i} - n b_{0} \sum_{i = 1}^{n} x_{i} - n b_{1} \sum_{i = 1}^{n} x^{2} = 0

$n\sum_{i=1}^n x_i \sum_{i=1}^ny_i -nb_0 \sum_{i=1}^n x_i-nb_1\sum_{i=1}^n x^2=0$

변형된 첫번째 수식에서 두번째 수식을 빼면 아래처럼 정리할 수 있습니다.

\sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i} - n \sum_{i = 1}^{n} x_{i} y_{i} - b_{1} (\sum_{i = 1}^{n} x_{i})^{2} + n b_{1} \sum_{i = 1}^{n} x^{2} = 0

$\sum_{i=1}^n x_i \sum_{i=1}^n y_i -n\sum_{i=1}^n x_iy_i-b_1(\sum_{i=1}^n x_i)^2 + nb_1\sum_{i=1}^n x^2 =0$

좌우변을 정리하면 아래와 같습니다.

b_{1} [n \sum_{i = 1}^{n} x_{i}^{2} - (\sum_{i = 1}^{n} x_{i})^{2}] = n \sum_{i = 1}^{n} x_{i} y_{i} - \sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i}

$b_1 \left[ n\sum_{i=1}^n x_i^2-(\sum_{i=1}^n x_i)^2 \right] = n\sum_{i=1}^n x_i y_i \ - \sum_{i=1}^n x_i \sum _{i=1}^ny_i$

이를 다시 정리하면 아래와 같습니다.

b_{1} = \frac{n \sum_{i = 1}^{n} x_{i} y_{i} - \sum_{i = 1}^{n} x_{i}^{[6]} \sum_{i = 1}^{n} y_{i}^{[7]}}{n \sum_{i = 1}^{n} x_{i}^{2} - (\sum_{i = 1}^{n} x_{i})^{2 [9]}}

$b_1 = \frac{n\sum_{i=1}^n x_i y_i \ - \sum_{i=1}^n x_i^{[6]} \sum_{i=1}^n y_i^{[7]}}{n\sum_{i=1}^n x_i^2-(\sum_{i=1}^n x_i)^{2[9]}}$

이렇게 변형된 수식으로부터 다음 과정이 유도됩니다.

과정 5 수식에서 분자 중 $\sum_{i=1}^n x_i$ [6로 표시한 부분]은 아래처럼 다시 쓸수 있습니다.

\sum_{i = 1}^{n} x_{i} = n \frac{1}{n} \sum_{i = 1}^{n} x_{i} = n \bar{x}

$\sum_{i=1}^n x_i =n\frac{1}{n}\sum_{i=1}^n x_i =n\bar{x} \$

$\sum_{i=1}^n y_i$ [7로 표시한 부분]도 아래처럼 다시 쓸 수 있습니다[7로 표시한 부분].

\sum_{i = 1}^{n} y_{i} = n \frac{1}{n} \sum_{i = 1}^{n} y_{i} = n \bar{y}

$\sum_{i=1}^n y_i =n\frac{1}{n}\sum_{i=1}^n y_i =n\bar{y} \$

따라서 과정 5의 수식에서 분자 부분만 다시 정리하면 최종적으로는 아래와 같이 쓸 수 있습니다.

n \sum_{i = 1}^{n} x_{i} y_{i} - \sum_{i = 1}^{n} x_{i} \sum_{i = 1}^{n} y_{i} = n \sum_{i = 1}^{n} x_{i} y_{i} - n^{2} \bar{x} \bar{y}

$n\sum_{i=1}^n x_i y_i \ - \sum_{i=1}^n x_i \sum_{i=1}^n y_i = n\sum_{i=1}^n x_i y_i - n^2\bar{x}\bar{y}$

분모 중에서 $(\sum_{i=1}^n x_i)^2$ [9로 표시한 부분]도 아래처럼 바꿔 쓸 수 있습니다.

\sum_{i = 1}^{n} x_{i} = n \frac{1}{n} \sum_{i = 1}^{n} x_{i} = n \bar{x} ∴ (\sum_{i = 1}^{n} x_{i})^{2} = n^{2} {\bar{x}}^{2}

$\sum_{i=1}^n x_i =n\frac{1}{n}\sum_{i=1}^n x_i =n\bar{x} \\ \therefore (\sum_{i=1}^n x_i)^2 = n^2\bar{x}^2$

과정 7과 9 결과를 이용해서 분모를 다시 정리하면 아래와 같습니다.

n \sum_{i = 1}^{n} x_{i}^{2} - (\sum_{i = 1}^{n} x_{i})^{2} = n \sum_{i = 1}^{n} x_{i}^{2} - n^{2} {\bar{x}}^{2}

$n\sum_{i=1}^n x_i^2-(\sum_{i=1}^n x_i)^2 = n\sum_{i=1}^n x_i^2 - n^2\bar{x}^2$

과정 8, 10의 결과를 이용해서 과정 5의 수식을 다시 써 보면 아래와 같습니다.

b_{1} = \frac{n \sum_{i = 1}^{n} x_{i} y_{i} - n^{2} \bar{x} \bar{y}}{n \sum_{i = 1}^{n} x_{i}^{2} - n^{2} {\bar{x}}^{2}}

$b_1 = \frac{ n\sum_{i=1}^n x_i y_i -n^2\bar{x}\bar{y}}{n\sum_{i=1}^n x_i^2 - n^2\bar{x}^2}$

과정 11의 수식에는 분자와 분모에 공통적으로 n을 포함하므로 정리하면 아래와 같습니다.

b_{1} = \frac{\sum_{i = 1}^{n} x_{i} y_{i} - n \bar{x} \bar{y}}{\sum_{i = 1}^{n} x_{i}^{2} - n {\bar{x}}^{2}}

$b_1 = \frac{ \sum_{i=1}^n x_i y_i -n\bar{x}\bar{y}}{\sum_{i=1}^n x_i^2 - n\bar{x}^2}$

과정 12의 수식 중 분모 부분을 아래처럼 변형하고,

\sum_{i = 1}^{n} x_{i}^{2} - n {\bar{x}}^{2} = \sum_{i = 1}^{n} x_{i}^{2} - 2 n {\bar{x}}^{2 [14]} + n {\bar{x}}^{2 [15]}

$\sum_{i=1}^n x_i^2 - n\bar{x}^2 = \sum_{i=1}^n x_i^2 -2n\bar{x}^{2[14]} +n\bar{x}^{2[15]}$

가운데의 $2n\bar{x}^2$ [14로 표시한 부분]은 아래처럼 변형됩니다.

2 n {\bar{x}}^{2} = 2 \bar{x} \sum_{i = 1}^{n} x_{i} (∵ \sum_{i = 1}^{n} x_{i} = n \bar{x})

$2n\bar{x}^2 = 2\bar{x}\sum_{i=1}^n x_i \\ (\because \ \sum _{i=1}^nx_i= n\bar{x})$

$n\bar{x}^2$ [15로 표시한 부분]은 $\sum _{i=1}^n\bar{x}^2$ 으로 고쳐쓸수 있습니다.
과정 14, 15의 결과를 이용, 분모 부분을 다시 쓰명 아래와 같습니다.

\sum_{i = 1}^{n} x_{i}^{2} - 2 n {\bar{x}}^{2} + n {\bar{x}}^{2} = \sum_{i = 1}^{n} x_{i}^{2} - 2 \bar{x} \sum_{i = 1}^{n} x_{i} + \sum_{i = 1}^{n} {\bar{x}}^{2}

$\sum_{i=1}^n x_i^2 -2n\bar{x}^2 +n\bar{x}^2 = \sum_{i=1}^n x_i^2 -2\bar{x}\sum_{i=1}^nx_i +\sum_{i=1}^n \bar{x}^2$

따라서 $\sum_{i=1}^n(x_i - \bar{x})^2$ 으로 고쳐 쓸수 있습니다.
과정 12의 수식에서 분자는 이렇게 변형할 수 있습니다.

\sum_{i = 1}^{n} x_{i} y_{i} - n \bar{x} \bar{y} = \sum_{i = 1}^{n} x_{i} y_{i} - 2 n \bar{x} {\bar{y}}^{[19]} + n \bar{x} {\bar{y}}^{[21]}

$\sum_{i=1}^n x_i y_i -n\bar{x}\bar{y} = \sum_{i=1}^n x_i y_i -2n\bar{x}\bar{y}^{[19]} + n\bar{x}\bar{y}^{[21]}$

$2n\bar{x}\bar{y}$ [19로 표시한 부분]은 아래처럼 변형할수 있습니다.

2 n \bar{x} \bar{y} = n \bar{x} {\bar{y}}^{[20]} + n \bar{x} {\bar{y}}^{[20]}

$2n\bar{x}\bar{y} = n\bar{x}\bar{y}^{[20]} + n\bar{x}\bar{y}^{[20]}$

$n\bar{x} \bar{y}$ [20로 표시한 부분]은 각각 $\bar{y}\sum_{i=1}^n x_i$ , $\bar{x}\sum_{i=1}^ny_i$ 로 쓸 수 있습니다.
$n\bar{x}\bar{y}$ [21로 표시한 부분]은 $\sum_{i=1}^n\bar{x}\bar{y}$ 로 쓸 수 있습니다.
그러면 분자의 형태는 아래와 같아 집니다.

\sum_{i = 1}^{n} x_{i} y_{i} - \bar{y} \sum_{i = 1}^{n} x_{i} - \sum_{i = 1}^{n} x_{i} y_{i} + \sum_{i = 1}^{n} \bar{x} \bar{y}

$\sum_{i=1}^n x_i y_i -\bar{y}\sum_{i=1}^n x_i - \sum_{i=1}^n x_i y_i + \sum_{i=1}^n \bar{x}\bar{y}$

따라서 분자도 $\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})$ 의 형태로 고쳐쓸 수 있습니다.
과정 17과 23의 정리를 이용, 과정 12의 수식을 바꿔주면 아래와 같습니다.

b_{1} = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y})}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}}

$b_1 = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i - \bar{x})^2}$

$b_0$ 에 관한 수식을 유도하는 과정

\sum_{i = 1}^{n} (y_{i} - b_{0} - b_{1} x_{i}) = 0

$\sum_{i=1}^n(y_i-b_0-b_1x_i) = 0$

\sum_{i = 1}^{n} y_{i} = n b_{0} + b_{1} \sum_{i = 1}^{n} x_{i} = 0 n b_{0} = \sum_{i = 1}^{n} y_{i} - b_{1} \sum_{i = 1}^{n} x_{i} b_{0} = \frac{1}{n} \sum_{i = 1}^{n} y_{i} - b_{1} \frac{\sum_{i = 1}^{n} x_{i}}{n} ∴ b_{0} = \bar{y} - b_{1} \bar{x}

$\sum_{i=1}^n y_i = nb_0 + b_1\sum_{i=1}^n x_i =0 \\n b_0 = \sum_{i=1}^n y_i- b_1\sum_{i=1}^n x_i \\ b_0 =\frac{1}{n} \sum_{i=1}^n y_i-b_1\frac{\sum_{i=1}^n x_i}{n} \\ \therefore b_0 = \bar{y}-b_1\bar{x}$

다시 결론만 요약하면 이렇습니다.
절편을 추정하는데 사용하는 수식은 아래와 같습니다.

b_{0} = \bar{y} - b_{1} \bar{x}

$\ b_{0} = \bar{y}-b_{1}\bar{x}$

기울기를 추정하는데 사용하는 수식은 아래와 같습니다.

b_{1} = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y})}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}}

$b_1 = \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}$

그러면 위 두 수식으로 직접 기울기와 절편을 계산해보고 추정한 회귀 모형의 기울기, 절편과 동일한지 살펴보겠습니다.
먼저 기울기 계산 결과입니다.

b1 <- sum((x-mean(x))*(y-mean(y)))/sum((x-mean(x))^2)
b1

## [1] 0.5275445

아래가 앞서 추정한 회귀 모형의 기울기 입니다.

m$coefficients[[2]]

## [1] 0.5275445

두 결과가 일치합니다. 그러면 절편도 확인해보겠습니다.

b0 <- mean(y)-b1*mean(x)
b0

## [1] -0.4529822

m$coefficients[[1]]

## [1] -0.4529822

R이 계산한 결과와 직접 계산한 결과가 일치했습니다.

지금까지 선형 회귀 모형의 회귀 계수를 추정하기 위해 경사 강하법 과 최소 제곱법을 사용해 봤습니다.
두 방법으로 계산한 결과는 동일했지만 경사 강하법은 프로그래밍이 필요하고 learning rate를 조절하는 등 조율이 사용이 어려웠습니다.
반대로 최소 제곱법은 프로그래밍이 필요없고 계산이 쉬웠습니다.
이러한 이유로 R을 포함한 대부분의 분석소프트웨어는 단순 회귀 분석의 경우 최소 제곱법을 사용합니다.
다음 챕터에서는 이렇게 만들어진 선형 회귀 모형을 평가하는 방법들을 살펴보겠습니다.

4.1 기울기와 절편의 표준오차

R로 회귀 모형을 만들면 아래 처럼 회귀 계수뿐 아니라 회귀 계수와 관련된 여러 추가 정보들을 확인할 수 있습니다.

summary.lm(m)

## 
## Call:
## lm(formula = y ~ x, data = df_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9099 -0.2229  0.2377  0.4435  1.0273 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.452982   0.422806  -1.071    0.325    
## x            0.527545   0.008091  65.200 8.75e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9476 on 6 degrees of freedom
## Multiple R-squared:  0.9986, Adjusted R-squared:  0.9984 
## F-statistic:  4251 on 1 and 6 DF,  p-value: 8.754e-10

단순 선형 회귀에서 얻을 수 있는 회귀 계수는 두개 뿐으로 R은 이 두 회귀 계수의 표준오차 정보를 제공합니다.
R은 표준오차(Standard error)를 Std.Error로 나타냅니다.
그러면 먼저 기울기의 표준오차를 살펴보겠습니다.

기울기의 표준오차는 추정된 기울기를 판단할 때 사용합니다.
그 표준오차가 기울기 보다 크면 좋은 추정치라고 할 수 없습니다.
확인된 기울기의 표준오차는 약 0.008로 상대적으로 작고 또 기울기보다도 훨씬 작습니다.

기울기의 표준오차를 계산하는 수식은 아래와 같습니다.

s (b) = \sqrt{\frac{1}{n - 2} \frac{\sum_{i = 1}^{n} (y_{i} - \hat{y_{i}})^{2}}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}}}

$s(b) =\sqrt{\frac{1}{n-2} \frac{\sum_{i=1}^n(y_{i}-\hat{y_{i}})^2}{\sum_{i=1}^n(x_{i}-\bar{x})^2}}$ 위 수식으로 기울기의 표준오차를 계산하면 아래와 같습니다.

sb1 <- sqrt(1/(8-2)*sum((df_data$y-m$fitted.values)^2)/sum((df_data$x-mean(df_data$x))^2))
sb1

## [1] 0.008091144

계산 결과는 R이 계산한 결과와 일치합니다.

기울기의 t-value도 확인할 수 있습니다.
t test의 t-value를 의미하는데 여기서는 $b_1$ 추정치를 $b_1=0$ 인 경우와 비교하여 그 차이가 유의한지 검정하는데 사용됩니다.
t-value를 계산하는 수식은 아래와 같습니다.

t = \frac{b_{1} - 0}{S E}

$t= \frac{b_1-0}{SE}$
SE는 기울기의 표준오차로 바로 위에서 계산했습니다.
직접 t-value를 계산하면 아래와 같습니다.

m$coefficients[[2]]/sb1

## [1] 65.20024

R이 계산한 t-value와 동일합니다. 이어서 위에서 계산된 t-value가 통계적으로 유의한지 살펴보겠습니다.

2*(1-pt(m$coefficients[[2]]/sb1, 6))

## [1] 8.753882e-10

이 결과도 R이 계산한 결과와 동일합니다.
계산된 값이 0.05보다 작기 때문에 둘의 차이는 통계적으로 유의하다고 판단할 수 있습니다.

같은 방법으로 절편의 표준오차와 t-value도 계산할 수 있습니다.
절편의 표준오차를 계산하는 수식은 아래와 같습니다.

\sqrt{\frac{\sum_{i = 1}^{n} (y - \hat{y})^{2}}{n - 2}} \sqrt{\frac{1}{n} + \frac{(\bar{x})^{2}}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}}}

$\sqrt{\frac{\sum_{i=1}^n{(y-\hat{y})^2}}{n-2}}\sqrt{\frac{1}{n}+ \frac{(\bar{x})^2}{\sum_{i=1}^n({x_{i}-\bar{x})^2}}}$

위 수식으로 계산한 결과인 절편의 표준오차는

sb0 <- sqrt(sum((y-m$fitted.values)^2)/(6))*sqrt((1/8)+(mean(x)^2)/sum((x-mean(x))^2))
sb0

## [1] 0.4228058

이고 R이 계산한 값과 동일합니다.

t-value는 아래와 같이 계산하고 R이 계산한 값과 같음을 알 수 있습니다.

m$coefficients[[1]]/sb0

## [1] -1.071372

마지막으로 추정된 절편이 통계적으로 유의한지 살펴보면 R이 계산한 값과 같고

2*(1-pt(abs(m$coefficients[[1]]/sb0), 6))

## [1] 0.3251937

그 값이 0.05보다 크다는 사실을 알 수 있습니다.
이 뜻은 절편의 추정량은 -0.453 이지만 통계적으로는 0과 다르지 않다고 할 수 없다는 뜻이 됩니다.
하지만 통계적으로 유의하지 않다고 하더라도 계산에 절편을 사용하지 않는다는 뜻은 아닙니다.

4.2 잔차의 표준오차

Residual standard error로 줄여서 RSE로 부르는데 이 수치가 작으면 작을 수록 실측값과 예측값이 잘 맞아 떨어진다는 뜻입니다.
따라서 모형의 정확도를 평가할 때 사용할 수 있고 특히 같은 데이터로 만들어진 서로 다른 회귀 모형이 존재할 때 비교 우위를
판단하는데 사용할 수 있는 척도가 되므로 알아 두면 좋은 개념입니다.

잔차의 표준오차를 계산하는 수식은 아래와 같습니다.

R S E = \sqrt{\frac{\sum_{i = 1}^{n} (y_{i} - f (x_{i}))^{2}}{n - 2}}

$RSE= \sqrt{\frac{\sum_{i=1}^{n}(y_i-f(x_i))^2}{n-2}}$

우리 회귀 모형의 RSE는 아래와 같습니다.

sqrt(sum(m$residuals^2)/6)

## [1] 0.9476273

R이 계산한 값과 동일합니다.

summary(m)$sigma

## [1] 0.9476273

RSE의 쓰임새를 살펴보기 위해 종종 사용되는 방법으로 가중치(weight)를 적용(여기서는 $\frac{1}{x^2}$ )한 모형과 그렇지 않은 모형을
비교해겠습니다. 먼저 가중치를 적용한 모형을 만듭니다.

m2 <- lm(y ~ x, df_data, weights = 1/x^2)
summary.lm(m2)

## 
## Call:
## lm(formula = y ~ x, data = df_data, weights = 1/x^2)
## 
## Weighted Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.029795 -0.013240 -0.007629  0.007874  0.042626 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01032    0.02868    0.36    0.731    
## x            0.50222    0.01171   42.89 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02624 on 6 degrees of freedom
## Multiple R-squared:  0.9967, Adjusted R-squared:  0.9962 
## F-statistic:  1840 on 1 and 6 DF,  p-value: 1.075e-08

가중치 적용으로 기울기와 절편이 모두 달라졌는데 주목할 만한 점은 RSE가 줄었다는 점입니다.
따라서 가중치를 적용한 모형이 그렇지 않은 모형보다 상대적으로 더 정확한 또는 예측을 잘하는 모형이라는 뜻이 됩니다.
두 모형을 한 그림으로 시각화 하면 아래와 같습니다.

ggplot(df_data, aes(x, y))+
  geom_point()+
  stat_function(fun = function(x, y)(y = m$coefficients[2]*x-m$coefficients[1]), col="red")+
  stat_function(fun = function(x, y)(y = m2$coefficients[2]*x-m2$coefficients[1]), col="blue")+
  theme_bw()+
  xlab("ng/mL") +
  ylab("Peack Response")+
  annotate("text", x = 100, y = 60, label = "Unweighted",parse = TRUE, col = "red") +
  annotate("text", x = 100, y = 45, label = "Weighted", parse = TRUE, col ="blue")

참고로 이렇게 가중치를 적용하는 경우도 회귀 계수를 추정하는데는 최소 제곱법을 사용합니다.
가중치를 사용하므로 수식은 아래 처럼 약간 변형됩니다.

b_{1} = \frac{\sum_{i = 1}^{n} w_{i} (x_{i} - {\bar{x}}_{w}) (y_{i} - {\bar{y}}_{w})}{\sum_{i = 1}^{n} w_{i} (x_{i} - {\bar{x}}_{w})^{2}}

$b_1 = \frac{\sum_{i=1}^{n}w_i(x_i-\bar{x}_w)(y_i-\bar{y}_w)}{\sum_{i=1}^{n}w_i(x_i-\bar{x}_w)^2}$ 위 수식에서 가중치는 아래처럼 표현하고,

w_{i} = \frac{1}{x_{i}^{2}}

$w_i = \frac{1}{x_i^2}$

가중치가 적용된 예측 변수의 평균은 아래 수식으로 계산하고,

{\bar{x}}_{w} = \frac{\sum_{i = 1}^{n} w_{i} x_{i}}{\sum_{i = 1}^{n} w_{i}}

$\bar{x}_w = \frac{\sum_{i=1}^{n}w_ix_i}{\sum_{i=1}^{n}w_i}$

가중치가 적용된 반응 변수(assay signal)의 평균도 같은 방식으로 계산합니다.

{\bar{y}}_{w} = \frac{\sum_{i = 1}^{n} w_{i} y_{i}}{\sum_{i = 1}^{n} w_{i}}

$\bar{y}_w = \frac{\sum_{i=1}^{n}w_iy_i}{\sum_{i=1}^{n}w_i}$

위 수식을 이용해서 직접 기울기의 회귀 계수를 계산해보겠습니다.

w <- 1/x^2
w.xbar <- sum(w*x)/sum(w)
w.ybar <- sum(w*y)/sum(w)
b1 <- sum(w*(x-w.xbar)*(y-w.ybar))/sum(w*(x-w.xbar)^2)
b1

## [1] 0.5022161

아래 R이 계산한 수치와 동일합니다.

m2$coefficients[2]

##         x 
## 0.5022161

RSE는 어떨까요?
가중치인 $\frac{1}{x^2}$ 을 곱하는 과정이 추가됐을 뿐 동일합니다.
아래 수식을 이용해서 직접 계산해보면 아래처럼 나오는데,

R S E = \sqrt{\frac{\sum_{i = 1}^{n} w_{i} (y_{i} - f (x_{i}))^{2}}{n - 2}}

$RSE= \sqrt{\frac{\sum_{i=1}^{n}w_i(y_i-f(x_i))^2}{n-2}}$

sqrt(sum(w*(y-m2$fitted.values)^2)/6)

## [1] 0.02624442

R이 계산한 결과와 역시 동일했습니다.

summary(m2)$sigma

## [1] 0.02624442

가중치가 적용된 경우, ordinary least squares method와 구분하기 위해서 weighted least squares method,
또는 줄여서 WLS라고 부르는데 원리는 같다는 것을 알 수 있었습니다.

결국 RSE가 더 작아졌다는 뜻은 그 만큼 더 정확해졌다는 뜻입니다.
따라서 보정하지 않은 모형보다 보정한 모형이 더 정확한 모형이라고 할 수 있습니다.
그런데 여기서 아무 근거 없이 $\frac{1}{x^2}$ 를 가중치로 사용한 것은 아닙니다.
예시로 사용된 데이터는 표준화 잔차가 x의 크기에 따라 커졌기 때문에 사용할 수 있었을 뿐입니다.
반대로 잔차가 x크기에 비례해서 커지는 경우 $\frac{1}{x^2}$ 를 가중치로 사용하면 모형의 정확도를 증가시킬 수 있다는 뜻이 됩니다.

아래에서 보면 가중치를 적용하지 않은 모형은 농도가 증가할 수록 표준화 잔차가 더 커졌습니다.

plot(m$fitted.values, rstandard(m))
abline(0,0)

반대로 가중치를 적용한 모형은 0을 중심으로 고른 표준화잔차 분포를 보입니다.

plot(m2$fitted.values, rstandard(m2))
abline(0,0)

따라서 주어진 데이터에 따라 적절한 가중치를 적용할 때만 모형의 정확도를 개선할 수 있다는 점이 중요한 점이라고 할 수 있겠습니다.

4.3 회귀 모형의 적합도

R은 피어슨 상관계수이고 이를 제곱한 값을 결정계수라고 부르는데 회귀 모형의 적합도(goodness of fit)를 판단하는데 사용합니다.
결정계수를 계산하는 수식은 아래와 같습니다.

S S T = \sum_{i = 1}^{n} (y_{i} - \bar{y})^{2}

$SST = \sum_{i=1}^{n}(y_i-\bar{y})^2$ SST(Total Sum of Squares)는 기본 모형이 설명하는 변동의 양으로 여기서 기본 모형은 기울기가 0인 모형입니다.
따라서 Y의 평균이 기본 모형이 됩입니다.

Y의 평균은 16.3625로 이때 각 실측값에서부터 평균까지 수직 선을 그리면 그 길이가 이 모형의 오차가 되고,
각 오차들을 제곱해서 합한 결과가 바로 SST입니다.

WLS 모형의 SST는 다음과 같이 계산합니다.

SST <- sum((w*(y-w.ybar)^2))
SST

## [1] 1.271126

ggplot(df_data, aes(x, y)) +
  geom_point(size=2, alpha=.4,
             position = position_jitter())+
  stat_function(fun = function(y) mean(df_data$y))+
  geom_segment(aes(xend=x, yend=mean(y)), lty=2)+
  theme_bw()

SSR(sum of squares regression)은 기본 모형과 회귀 모형의 간 차이의 제곱합을 의미합니다.

S S R = \sum_{i = 1}^{n} (\hat{y_{i}} - \bar{y})^{2}

$SSR = \sum_{i=1}^{n}(\hat{y_i}-\bar{y})^2$

WLS 모형의 SSR은 아래처럼 계산합니다.

SSR <- sum(w*(m2$fitted.values-w.ybar)^2)
SSR

## [1] 1.266993

ggplot(df_data, aes(x, y)) +
  geom_point(size=2, alpha=.4,
             position = position_jitter())+
  stat_function(fun = function(x, y)(y = m2$coefficients[2]*x-m2$coefficients[1]), col="blue")+
  stat_function(fun = function(y) mean(df_data$y)) +
  geom_segment(aes(x = x, 
                   xend = x, 
                   y = mean(y), 
                   yend = predict(m2, newdata=df_data)), lty=2) +
  theme_bw()

SSE(Sum of squares error)는 이미 알고 있지만 다시 설명하면 예측값과 실측값의 차이의 제곱합을 의미합니다.

S S E = \sum_{i = 1}^{n} (\hat{y_{i}} - y_{i})^{2}

$SSE = \sum_{i=1}^{n}(\hat{y_i}-y_i)^2$

WLS 모형에서 SSE는 아래와 같이 계산합니다.

SSE <- sum(w*(m2$fitted.values-y)^2)
SSE

## [1] 0.004132618

ggplot(df_data, aes(x, y)) +
  geom_point(size=2, alpha=.4,
             position = position_jitter())+
  stat_function(fun = function(x, y)(y = m2$coefficients[2]*x-m2$coefficients[1]), col="blue")+
  geom_segment(aes(xend=x, yend=m2$fitted.values), lty=2)+
  theme_bw()

SST는 SSR과 SSE를 더한 것과 동일한데 아래와 같이 증명됩니다.

SST == SSE+SSR

## [1] TRUE

결정계수를 계산하는 방법은 아래와 같습니다.

R^{2} = \frac{S S R}{S S T}

$R^{2} = \frac{SSR}{SST}$

SSR/SST

## [1] 0.9967489

summary(m2)$r.squar

## [1] 0.9967489

R이 계산한 결과와 직접 계사한한 결과가 같았습니다.

또는 SST=SSR+SSE를 이용해서 아래처럼 변형해서 계산할 수 있습니다.

R^{2} = \frac{S S T - S S E}{S S T}

$R^{2} = \frac{SST-SSE}{SST}$

R^{2} = 1 - \frac{S S E}{S S T}

$R^{2} = 1- \frac{SSE}{SST}$

1-SSE/SST

## [1] 0.9967489

summary(m2)$r.squar

## [1] 0.9967489

그런데 결정계수는 정량 분석에 사용하기 좋은 척도는 아닙니다.
회귀 모형으로 실측값(assay signal)을 농도로 바꿔보면 아래와 같습니다.

back.x <- round((y-m$coefficients[[1]])/m$coefficients[[2]],3)
back.x2 <- round((y-m2$coefficients[[1]])/m2$coefficients[[2]],3)

가중치를 적용하지 않은 모형의 back calculation 농도는
1.806, 2.944, 4.46, 8.441, 15.644, 31.377, 60.38, 129.947
로 낮은 농도에서 정확성이 떨어집니다.

반대로 가중치를 적용한 모형의 back calculation 농도는
0.975, 2.17, 3.763, 7.944, 15.511, 32.037, 62.502, 135.578
로 낮은 농도에서도 상대적으로 정확성이 높습니다.

더 비교하기 쉽도록 Accuracy(%) 단위로 바꿔서 표현해보면, 가중치를 적용하지 않은 모형:

round((back.x-x)/x*100, 1)

## [1] 80.6 47.2 11.5  5.5 -2.2 -1.9 -5.7  1.5

가중치를 적용한 모형:

round((back.x2-x)/x*100, 1)

## [1] -2.5  8.5 -5.9 -0.7 -3.1  0.1 -2.3  5.9

즉 가중치를 적용하지 않은 모형은 낮은 농도의 정확도가 매우 부정확하다는 사실을 알 수 있습니다.
그러나 두 모형의 결정계수는 모두 0.99 이상입니다.

가중치를 적용하지 않은 모형의 결정계수:

summary(m)$r.squared

## [1] 0.9985906

가중치를 적용한 모형의 결정계수:

summary(m2)$r.squared

## [1] 0.9967489

오히려 가중치를 적용한 모형의 결정계수가 조금 더 작습니다.

결정계수는 예측변수의 숫자가 늘어나면 커지는 성질이 있습니다.
이를 보완하기 위해 Adjusted R squared를 사용하는게 더 일반적입니다.
Adjusted R squared 계산 방법은 아래와 같습니다.

A d j u s t e d R^{2} = 1 - \frac{\frac{S S E}{n - k - 1}}{\frac{S S T}{n - 1}}

$Adjusted \ R^2 \ = 1-\frac{\frac{SSE}{n-k-1}}{\frac{SST}{n-1}}$

K는 예측변수의 수인데 농도라는 단일 예측변수를 사용했으므로 k=1입니다.
n은 표본의 수로 8개의 농도를 측정했기 때문에 n=8입니다.
아래 직접 계산한 결과와 R이 계산한 결과를 비교하면 역시 동일하다는 것을 알 수 있습니다.

1-(SSE/(8-1-1))/(SST/(8-1))

## [1] 0.996207

summary(m2)$adj.r.squared

## [1] 0.996207

또 아래 가중치를 적용하지 않은 모형의 Adjusted R squared 수치와 비교하면 역시 조금 작아졌습니다.

summary(m)$adj.r.squared

## [1] 0.9983557

하지만 설명한바와 같이 결정계수 자체의 크기는 이미 1에 가까운 수치이고 가중치를 적용한 모형이 가중치를 적용하지 않은 모형보다
결정계수 크기가 더 작습니다.
따라서 결정계수의 수치만으로 정량 분석에 사용할 회귀 모형의 적합도를 판단하는 것은 적절하지 않다는 사실을 알 수 있습니다.

4.4. F-statistic

마지막입니다.

summary.lm(m2)

## 
## Call:
## lm(formula = y ~ x, data = df_data, weights = 1/x^2)
## 
## Weighted Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.029795 -0.013240 -0.007629  0.007874  0.042626 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01032    0.02868    0.36    0.731    
## x            0.50222    0.01171   42.89 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02624 on 6 degrees of freedom
## Multiple R-squared:  0.9967, Adjusted R-squared:  0.9962 
## F-statistic:  1840 on 1 and 6 DF,  p-value: 1.075e-08

선형 회귀 모명을 출력하면 마지막에 F-statistic이 보이는데 이는 F비(ratio)로 아래 처럼 계산합니다.

F = \frac{M S R}{M S E}

$F \ = \frac{MSR}{MSE}$

조금 더 살펴보면,
MSR은 SSR를 예측변수의 수로 나눈 값입니다. 그래서 Mean square regression이라고 부릅니다.
SSR은 앞에서도 살펴봤지만 예측값과 평균 assay signal 간 차이의 제곱합입니다.
또 예측변수로는 농도만 사용되었으므로 k=1이라 이 경우 MSR은 SSR이 됩니다.

M S R = \frac{S S R}{k}

$MSR = \frac{SSR}{k}$

MSR <- SSR/1
MSR

## [1] 1.266993

MSE는 Mean squared error로 SSE를 표본 수에서 예측 변수의 수와 1을 뺀 값으로 나눈 값입니다.
따라서 SSE를 6(8-1-1 = 6)으로 나눈 값 입니다.
SSE는 설명한대로 예측값과 실측값 간 차이의 제곱합이였습니다.
MSE는 사실 cost function을 설명하면서 잠깐 언급되기도 했습니다.

M S E = \frac{S S E}{n - k - 1}

$MSE = \frac{SSE}{n-k-1}$

MSE <- SSE/6
MSE

## [1] 0.0006887697

MSR과 MSE를 이용해서 F를 계산하면 아래와 같이 R이 계산한 결과와 동일한 것을 알 수 있습니다.

MSR/MSE

## [1] 1839.502

F비의 의미는 따라서 MSR과 MSE의 비율로 모형이 설명하는 변동과 모형이 설명하지 않는 변동의 비율인 셈입니다.
따라서 F비가 클 수록 모형이 데이터를 더 잘 설명한다고 해석합니다.
이 F비도 통계적으로 유의한지 검정이 가능합니다.

qf(1 - 0.05, 1, 6)

## [1] 5.987378

유의수준 0.05, 각 자유도 1, 6인 경우 F비(임계값)는 약 5.99으로 우리 모형의 F비 통계량 1840 보다 훨씬 작습니다.
따라서 우리 모형의 F비는 통계적으로 유의하다고 판단합니다.

t-value는 아래처럼 계산하는데 R이 계산한 결과와 동일합니다.
마찬가지로 0.05보다 매우 작은 수치이므로 F비는 통계적으로 유의하다고 결론낼 수 있습니다.

pf(MSR/MSE, 1, 6, lower.tail = FALSE)

## [1] 1.075202e-08

지금까지 선형 회귀 분석과 그 방법을 공유하는 linear regression curve fitting을 살펴봤습니다.
서로 이름만 달랐을 뿐 같은 개념임을 알 수 있었고 경사강하법과 최소 제곱법을 비교해보기도 했습니다.
그러나 오늘 다룬 내용이 선형 회귀의 모든 부분은 아닙니다.
가령 선형대수를 이용하는 방법은 더 많은 수학적 배경지식이 필요해서 다루지 않았습니다.
그만큼 분석 업무는 많은 부분을 수학과 통계학에 의존하고 있다고 해도 과언이 아닐 것입니다.
다음 시간에는 Ligand binding assay, cell based assay에 사용되는 Non-linear regression을 다뤄보겠습니다.